Packages
`%nin%` = Negate(`%in%`)
if (!require("tidyverse")) install.packages("tidyverse")
library("tidyverse")# for dplyr workflow of calculations
if (!require("UsingR")) install.packages("UsingR")
library("UsingR")
if (!require("lme4")) install.packages("lme4")
library("lme4")
if (!require("lmerTest")) install.packages("lmerTest")
library("lmerTest")
if (!require("broom")) install.packages("broom")
library("broom") # lmer model export
if (!require("broom.mixed")) install.packages("broom.mixed")
library("broom.mixed") # lmer model export
if (!require("ggplot2")) install.packages("ggplot2")
library("ggplot2")
if (!require("viridis")) install.packages("viridis")
library("viridis")
if (!require("RColorBrewer")) install.packages("RColorBrewer")
library("RColorBrewer")
if (!require("ggpubr")) install.packages("ggpubr")
library("ggpubr") # for multi-ggplot figs
if (!require("Hmisc")) install.packages("Hmisc")
library("Hmisc") # correlation matrix
if (!require("emmeans")) install.packages("emmeans")
library("emmeans") # marginal means & confintsSet Colors & Themes
savs_ggplot_theme <- theme_bw() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text.align = 0,
legend.position = "bottom"
#plot.margin = margin(t = 6, r = 6, b = 0, l = 1, unit = "pt")
)
my_vir_cols_7 <- viridis(n = 7,
#alpha,
#begin, # determine how dark it can go
end = 0.8, # makes it so it cannot get too light in color
#direction = -1,
option = "mako"
)
my_vir_cols_41 <- viridis(n = 4, # how many different colors do I need?
alpha = 1, # transparency
begin = 0, # how dark it's allowed to get
end = 0.95, # how light it's allowed to get
direction = -1
)
my_vir_cols_42 <- magma(n = 4, # how many different colors do I need?
alpha = 1, # transparency
begin = 0, # how dark it's allowed to get
end = 0.95, # how light it's allowed to get
direction = -1
)
my_vir_cols_43 <- plasma(n = 4, # how many different colors do I need?
alpha = 1, # transparency
begin = 0, # how dark it's allowed to get
end = 0.95, # how light it's allowed to get
direction = -1
)
my_vir_cols_44 <- cividis(n = 4, # how many different colors do I need?
alpha = 1, # transparency
begin = 0, # how dark it's allowed to get
end = 0.95, # how light it's allowed to get
direction = -1
)
my_vir_cols_45 <- mako(n = 4, # how many different colors do I need?
alpha = 1, # transparency
begin = 0, # how dark it's allowed to get
end = 0.95, # how light it's allowed to get
direction = -1
)
my_brew_cols_4 <- brewer.pal(9, "YlGnBu")[c(3,5,7,9)]Background & Introduction
Telemetry response variables to test vs hydric physiology - Microhabitat use = proportion of observations doing x (above, partial shade, full shade, burrow) - Time of onset of estivation (once you no longer ever get aboveground obs)
Load Data
Lizard Data
This is the same dataframe used in ‘analysis_1’.
apr_may_dates <- as.Date(c("2021-04-23", "2021-04-24", "2021-04-25"))
may_dates <- as.Date(c("2021-05-07", "2021-05-08", "2021-05-09"))
liz_dat <- read_rds("./data/G_sila_clean_full_dat.RDS") %>%
# we only use data for radio collared lizards in this analysis
dplyr::filter(complete.cases(radio_collar_serial)) %>%
# add date chunks to help look at values over time
mutate(date_chunks = factor(case_when(capture_date %in% apr_may_dates ~ "Apr 26 - May 6",
capture_date %in% may_dates ~ "May 10 - 20")),
radio_collar_serial = factor(radio_collar_serial))
summary(liz_dat)## capture_date_time capture_date radio_collar_serial
## Min. :2021-04-23 10:39:00.00 Min. :2021-04-23 229-069: 3
## 1st Qu.:2021-04-23 14:54:30.00 1st Qu.:2021-04-23 245-742: 3
## Median :2021-04-24 10:50:30.00 Median :2021-04-25 245-744: 3
## Mean :2021-04-28 14:03:08.27 Mean :2021-05-12 245-745: 3
## 3rd Qu.:2021-05-07 11:41:30.00 3rd Qu.:2021-05-08 245-748: 3
## Max. :2021-05-08 14:19:00.00 Max. :2021-07-14 252-883: 3
## NA's :18 (Other):58
## PIT_tag_ID individual_ID sex_M_F mass_g SVL_mm
## Length:76 F-12 : 3 Female:32 Min. :25.00 Min. : 85.0
## Class :character M-09 : 3 Male :44 1st Qu.:33.70 1st Qu.: 96.0
## Mode :character M-10 : 3 Median :38.00 Median :100.0
## M-11 : 3 Mean :39.20 Mean :101.7
## M-19 : 3 3rd Qu.:43.45 3rd Qu.:108.0
## M-20 : 3 Max. :58.80 Max. :122.0
## (Other):58
## tmt tmt_mass_change_g capture_to_msmt hematocrit_percent
## Water Tmt:38 Min. :-3.0000 Min. : 11.0 Min. :17.00
## Sham Tmt :38 1st Qu.:-1.0000 1st Qu.: 62.5 1st Qu.:30.00
## No Tmt : 0 Median : 0.0000 Median : 152.0 Median :34.00
## Mean :-0.2262 Mean : 401.1 Mean :34.35
## 3rd Qu.: 0.0000 3rd Qu.: 265.2 3rd Qu.:38.50
## Max. : 3.0000 Max. :1665.0 Max. :58.00
## NA's :15 NA's :18 NA's :1
## cloacal_temp_C CEWL_g_m2h msmt_temp_C msmt_RH_percent
## Min. :24.00 Min. : 0.650 Min. :19.07 Min. :12.20
## 1st Qu.:30.00 1st Qu.: 7.989 1st Qu.:28.02 1st Qu.:14.79
## Median :32.00 Median : 9.620 Median :29.31 Median :20.38
## Mean :31.92 Mean : 9.854 Mean :29.18 Mean :22.70
## 3rd Qu.:34.00 3rd Qu.:12.324 3rd Qu.:31.74 3rd Qu.:27.49
## Max. :36.00 Max. :21.673 Max. :33.55 Max. :40.22
##
## e_s_kPa msmt_VPD_kPa osmolality_mmol_kg_mean month week
## Min. :2.205 Min. :1.415 Min. :309.3 April:39 1 :39
## 1st Qu.:3.784 1st Qu.:2.487 1st Qu.:345.0 May :23 3 :23
## Median :4.076 Median :3.229 Median :361.0 July :14 12:14
## Mean :4.100 Mean :3.212 Mean :363.6
## 3rd Qu.:4.685 3rd Qu.:3.982 3rd Qu.:377.9
## Max. :5.188 Max. :4.319 Max. :436.5
## NA's :1
## week_num ultrasound_date gravid_Y_N number_eggs
## Min. : 1.000 Min. :2021-04-24 Gravid :16 Min. :3.000
## 1st Qu.: 1.000 1st Qu.:2021-04-24 Not Gravid:12 1st Qu.:3.750
## Median : 1.000 Median :2021-04-24 NA's :48 Median :4.000
## Mean : 3.632 Mean :2021-04-29 Mean :3.812
## 3rd Qu.: 3.000 3rd Qu.:2021-05-08 3rd Qu.:4.000
## Max. :12.000 Max. :2021-05-08 Max. :5.000
## NA's :48 NA's :60
## largest_egg_diameter_cm largest_egg_circumference_cm stage
## Min. :0.5700 Min. :1.710 small round: 1
## 1st Qu.:0.8425 1st Qu.:2.385 large round: 9
## Median :0.9950 Median :2.885 soft : 2
## Mean :1.1131 Mean :3.148 soft oblong: 1
## 3rd Qu.:1.4200 3rd Qu.:3.985 firm oblong: 2
## Max. :1.8000 Max. :4.740 hard oblong: 1
## NA's :60 NA's :60 NA's :60
## dev_point SMI date_chunks
## Min. :1.00 Min. :23.67 Apr 26 - May 6:39
## 1st Qu.:2.00 1st Qu.:32.29 May 10 - 20 :23
## Median :2.00 Median :34.40 NA's :14
## Mean :2.75 Mean :34.77
## 3rd Qu.:3.25 3rd Qu.:37.43
## Max. :5.00 Max. :45.63
## NA's :60liz_dat_apr_may_only <- liz_dat %>%
dplyr::filter(complete.cases(date_chunks))
liz_dat_july_only <- liz_dat %>%
dplyr::filter(month == "July")
liz_IDs_only <- liz_dat %>%
group_by(radio_collar_serial, individual_ID, sex_M_F) %>%
summarise(n = n()) %>%
dplyr::select(-n)## `summarise()` has grouped output by 'radio_collar_serial', 'individual_ID'. You
## can override using the `.groups` argument.Tsurf + Wrangling
We measured the surface body temperature of lizards throughout the study using Holohil temperature-sensitive transmitters, which is the data we load in this section.
variables in this dataframe:
- DateTime = date and time the temperature was logged
- Serial = the radio-collar serial number for that observation (to be used to link to specific individuals)
- Freq = what radio frequency the collar was on
- pp = the pulse period recorded by the radio-collar receiver (pulse period is a proxy for temperature)
- ppTempC = the pp value converted into temperature in degrees Celsius using third order polynomial functions estimated from the reference transformations provided by Holohil
temps <- read.csv("./data/telemetry/tower_temps_all.csv") %>%
# properly format data classes
mutate(date_only = as.Date(substr(DateTime, 1, 10), format = "%Y-%m-%d"),
hour_of_day = as.numeric(substr(DateTime, 12, 13)),
DateTime = as.POSIXct(DateTime, format = "%Y-%m-%d %H:%M:%S"),
Serial = factor(paste(substr(Serial, 1, 3),
substr(Serial, 4, 7), sep = "-")),
Freq = round(Freq, 2))
# save freq-serial link
freq_serial <- temps %>%
group_by(Freq, Serial) %>%
summarise(n = n()) #%>%## `summarise()` has grouped output by 'Freq'. You can override using the
## `.groups` argument.# check that there is 1 serial per frequency
# group_by(Freq) %>%
# summarise(n = n())
# look at temps dataframe
summary(temps)## DateTime Serial Freq
## Min. :2021-04-25 13:11:35.00 245-762: 5952 Min. :150.0
## 1st Qu.:2021-05-02 18:18:19.75 245-756: 5610 1st Qu.:150.4
## Median :2021-05-11 07:58:47.00 245-753: 5606 Median :151.1
## Mean :2021-05-15 02:34:25.98 229-078: 5481 Mean :151.0
## 3rd Qu.:2021-05-17 22:48:35.50 252-886: 4899 3rd Qu.:151.5
## Max. :2021-07-13 13:11:42.00 229-069: 4728 Max. :152.0
## (Other):75176
## ppTempC pp date_only hour_of_day
## Min. :-23.99 Min. :1001 Min. :2021-04-25 Min. : 0.00
## 1st Qu.: 31.80 1st Qu.:1565 1st Qu.:2021-05-02 1st Qu.:10.00
## Median : 36.26 Median :1703 Median :2021-05-11 Median :13.00
## Mean : 35.20 Mean :1734 Mean :2021-05-14 Mean :12.69
## 3rd Qu.: 39.06 3rd Qu.:1885 3rd Qu.:2021-05-17 3rd Qu.:16.00
## Max. : 76.93 Max. :2799 Max. :2021-07-13 Max. :23.00
## unique(temps$Serial)## [1] 229-059 229-060 229-063 229-065 229-066 229-069 229-070 229-072 229-073
## [10] 229-074 229-077 229-078 229-079 229-080 245-741 245-742 245-745 245-748
## [19] 245-751 245-752 245-753 245-754 245-756 245-757 245-759 245-760 245-761
## [28] 245-762 252-881 252-882 252-883 252-884 252-885 252-886 252-887
## 35 Levels: 229-059 229-060 229-063 229-065 229-066 229-069 229-070 ... 252-887Subset the temperature data so that we do not use the temperatures from during up to 24h after the capturing and releasing lizards: NOTE I will also need to remove any data for lizards that went missing, and for when we re-assigned collars. These two radio collar serial numbers will be problematic. I need to go in to the Radio-collar temp data and add the “a” for data after the date we switched the collar to a new lizard.
to do - 245-762-a - 229-070-a
Also add columns for time chunks to make averaging easier later.
# capture weekend dates
cap_dates <- as.Date(c("2021-04-23", "2021-04-24", "2021-04-25", # April
"2021-05-07", "2021-05-08", "2021-05-09"), # May
format = "%Y-%m-%d")
# grouping notes by date
group_dates <- data.frame(date_only = c(seq(as.Date("2021-04-26"),
as.Date("2021-07-13"), by = 1))) %>%
dplyr::filter(date_only %nin% cap_dates) %>%
mutate(date_chunks = factor(c(rep("Apr 26 - May 6", 11),
rep("May 10 - 20", 11),
rep("May 21 - 31", 11),
rep("Jun 1 - 11", 11),
rep("Jun 12 - 22", 11),
rep("Jun 23 - Jul 3", 11),
rep("Jul 4 - 13", 10)),
levels = c("Apr 26 - May 6", "May 10 - 20",
"May 21 - 31", "Jun 1 - 11", "Jun 12 - 22",
"Jun 23 - Jul 3", "Jul 4 - 13")))
# subset dataset and add grouping notes
temp_sub <- temps %>%
dplyr::filter(date_only %nin% cap_dates) %>%
left_join(group_dates, by = 'date_only')
summary(temp_sub)## DateTime Serial Freq
## Min. :2021-04-26 00:01:18.0 245-762: 5559 Min. :150.0
## 1st Qu.:2021-05-02 12:53:38.0 245-756: 5195 1st Qu.:150.4
## Median :2021-05-12 07:38:43.0 245-753: 5183 Median :151.1
## Mean :2021-05-15 19:07:48.7 229-078: 4971 Mean :151.0
## 3rd Qu.:2021-05-22 11:00:08.0 252-886: 4487 3rd Qu.:151.5
## Max. :2021-07-13 13:11:42.0 229-069: 4457 Max. :152.0
## (Other):68431
## ppTempC pp date_only hour_of_day
## Min. :-23.99 Min. :1001 Min. :2021-04-26 Min. : 0.00
## 1st Qu.: 32.01 1st Qu.:1560 1st Qu.:2021-05-02 1st Qu.:10.00
## Median : 36.39 Median :1697 Median :2021-05-12 Median :13.00
## Mean : 35.33 Mean :1728 Mean :2021-05-15 Mean :12.61
## 3rd Qu.: 39.15 3rd Qu.:1877 3rd Qu.:2021-05-22 3rd Qu.:16.00
## Max. : 76.93 Max. :2798 Max. :2021-07-13 Max. :23.00
##
## date_chunks
## Apr 26 - May 6:39563
## May 10 - 20 :33856
## May 21 - 31 : 9545
## Jun 1 - 11 : 6263
## Jun 12 - 22 : 3431
## Jun 23 - Jul 3: 2700
## Jul 4 - 13 : 2925I also need to remove outliers, which are pretty standard erroneous points when using temperature-sensitive telemetry with devices external to the animal. I will follow Nicole Gaudenti’s 2021 paper, and exclude “any outliers greater than two standard deviations away from each lizard’s mean Tb”
, as these likely represented glitches in the data acquisition system; such outliers were uncommon (<5% of data points).
Look at data distribution:
boxplot(temp_sub$ppTempC)
hist(temp_sub$ppTempC)
ggplot(temp_sub) +
geom_point(aes(x = Serial,
y = ppTempC,
color = Serial)) +
theme_classic()
ggplot(temp_sub) +
geom_boxplot(aes(x = Serial,
y = ppTempC,
color = Serial)) +
theme_classic()
Get mean and SD for each lizard’s Tb data, then use it to filter out extreme values:
temp_clean <- temp_sub %>%
group_by(Serial) %>%
#summarise(outs = boxplot.stats(ppTempC)$out) # less outliers found this way
mutate(mean_Tb = mean(ppTempC),
sd_Tb = sd(ppTempC)) %>%
mutate(accept_low = mean_Tb - (2*sd_Tb),
accept_high = mean_Tb + (2*sd_Tb)) %>%
dplyr::filter(ppTempC < accept_high & ppTempC > accept_low) %>%
dplyr::select(-mean_Tb, -sd_Tb, -accept_low, -accept_high)
summary(temp_clean)## DateTime Serial Freq
## Min. :2021-04-26 00:01:18.00 245-762: 5337 Min. :150.0
## 1st Qu.:2021-05-02 16:03:25.00 245-756: 4984 1st Qu.:150.4
## Median :2021-05-12 10:31:25.00 245-753: 4974 Median :151.1
## Mean :2021-05-16 01:47:49.74 229-078: 4686 Mean :151.0
## 3rd Qu.:2021-05-22 14:43:39.00 229-069: 4294 3rd Qu.:151.5
## Max. :2021-07-13 13:11:42.00 252-886: 4251 Max. :152.0
## (Other):65267
## ppTempC pp date_only hour_of_day
## Min. : 1.67 Min. :1109 Min. :2021-04-26 Min. : 0.00
## 1st Qu.:32.44 1st Qu.:1563 1st Qu.:2021-05-02 1st Qu.:10.00
## Median :36.47 Median :1693 Median :2021-05-12 Median :13.00
## Mean :35.56 Mean :1720 Mean :2021-05-15 Mean :12.74
## 3rd Qu.:39.09 3rd Qu.:1861 3rd Qu.:2021-05-22 3rd Qu.:16.00
## Max. :53.17 Max. :2798 Max. :2021-07-13 Max. :23.00
##
## date_chunks
## Apr 26 - May 6:36972
## May 10 - 20 :32632
## May 21 - 31 : 9293
## Jun 1 - 11 : 6098
## Jun 12 - 22 : 3316
## Jun 23 - Jul 3: 2638
## Jul 4 - 13 : 2844What proportion of points did we remove?
(nrow(temp_sub) - nrow(temp_clean))/nrow(temp_sub)## [1] 0.04568444.6%
Plot again:
ggplot(temp_clean) +
geom_point(aes(x = Serial,
y = ppTempC,
color = Serial)) +
theme_classic()
ggplot(temp_clean) +
geom_boxplot(aes(x = Serial,
y = ppTempC,
color = Serial)) +
theme_classic()
Individual IQRs definitely look reasonable.
Microhabitat Use
This data was collected by radio-tracking lizards, finding them, then recording what type of microhabitat they were in. We have to read-in a few files, because at first the formatting was not consistent, then put the dataframes all together.
variables in this dataframe:
- Serial = the radio-collar serial number for that observation (to be used to link to specific individuals)
- DateTime = date and time the temperature was logged
- microhabitat (was Micro_1) = which microhabitat lizards were observed in
- Breed_Col = breeding coloration of lizards (we will not use)
- bpmTempC = the temperature in Celsius from the lizards’ collar, based on the handheld receivers used to locate them, which record the pulse interval in beats per minute rather than the pulse period of the tower receiver for other Tb data
- dwsTempOut = the temperature recorded by a weather station at the study site; these values were interpolated to the specific time of the observation since the weather station was only recording data every two hours ModTempC = the ground surface temperature estimated by a temperature model (Ian Axsom) for the location and time that each observation occurred; this is an estimate of the operative temperature available to these lizards
This is the nice, clean data that was collected early May onwards:
microhab_clean <- read.csv("./data/telemetry/behavior_observations.csv") %>%
# force NAs on, then remove, observations for the "uncollared" lizard
mutate(Serial_no = as.numeric(substr(Serial, 1, 5))) %>%
dplyr::filter(complete.cases(Serial_no)) %>%
# properly format data classes
mutate(date_only = as.Date(substr(DateTime, 1, 10), format = "%Y-%m-%d"),
hour_of_day = as.numeric(substr(DateTime, 12, 13)),
DateTime = as.POSIXct(DateTime, format = "%Y-%m-%d %H:%M:%S"),
Serial = factor(paste(paste(substr(as.character(Serial), 1, 3),
substr(as.character(Serial), 4, 7), sep = "-"))),
microhabitat = factor(Micro_1,
labels = c("Burrow", # keep
"Open", # change from Burrow slide
"Open", # change from Burrow-Partial
"Burrow", # change from Burrow/ephedra
"Burrow", # change from Burrow/shrub
"", ### change from Dead/Missing ### ### ### ###
"Open", # change from Grass trims
"", # change from IDK
"", # change from No Visual
"Open", # keep
"Shade - Full", # keep
"Shade - Partial", # keep
"Shade - Partial", # remove notes
"Shade - Partial", # remove notes
"Shade - Partial", # remove notes
"Shade - Partial", # remove notes
"Shade - Partial", # remove notes
"Shade - Partial", # remove notes
"Shade - Partial" # remove notes
))
) %>%
# for some reason, it only works to fix the order separately
mutate(microhabitat = factor(microhabitat,
levels = c("Open",
"Shade - Partial",
"Shade - Full",
"Burrow"
#"No Visual", "Dead/Missing"
))) %>%
dplyr::select(Serial, DateTime, date_only, hour_of_day, microhabitat,
bpmTempC, dwsTempOut, ModTempC) %>%
# add date chunk grouping notes
dplyr::filter(date_only %nin% cap_dates)## Warning in mask$eval_all_mutate(quo): NAs introduced by coercionThis is the messy data collected in late April and very early May. Unfortunately, the following radio-frequencies were unable to be linked to a serial number: 150.14, 150.19, 150.21, and 150.28. But, that only applied to 11 observations.
microhab_early0 <- read.csv("./data/telemetry/behavior_observations_early_0.csv") %>%
dplyr::select(date_only = date, time,
microhabitat = micro,
Serial = serial,
Freq = freq, freq_notes = X)
# this one must be split by whether freq or serial number of radio collars recorded
microhab_early0_freq_only <- microhab_early0 %>%
dplyr::filter(!complete.cases(Serial)) %>%
dplyr::select(-Serial)
microhab_early0_serial <- microhab_early0 %>%
dplyr::filter(complete.cases(Serial)) %>%
mutate(Serial = factor(paste(substr(Serial, 1, 3),
substr(Serial, 4, 7), sep = "-")))
# load other csv, and add the dfs above to it
microhab_early_all <- read.csv("./data/telemetry/behavior_observations_early_1.csv") %>%
dplyr::select(date_only = date, time,
microhabitat = micro,
Freq = freq, freq_notes) %>%
rbind(microhab_early0_freq_only) %>%
mutate(Freq = round(Freq, 2)) %>%
left_join(freq_serial, by = 'Freq') %>%
dplyr::select(-n) %>%
rbind(microhab_early0_serial) %>%
mutate(DateTime = as.POSIXct(paste(date_only, time, sep = " "),
format = "%m/%d/%y %H:%M"),
hour_of_day = as.numeric(substr(as.character(DateTime), 12, 13)),
date_only = as.Date(date_only, format = "%m/%d/%y"),
microhabitat = factor(str_to_title(microhabitat))) %>%
dplyr::filter(complete.cases(Serial, microhabitat)) %>%
# some of these were "focal observations" taken every few minutes
# to be consistent with other obs, take first point of each focal
group_by(Serial, date_only) %>%
#summarise(min(date), max(date), min(time), max(time)) # check that data is one chunk
# might remove a few usable data, but OK
slice(1) %>%
# remove these to make it easy to bind to clean one
dplyr::select(-time, -Freq, -freq_notes) %>%
# add these to make it easy to bind to clean one
mutate(bpmTempC = NA, dwsTempOut = NA, ModTempC = NA)
summary(microhab_early_all) ## date_only microhabitat Serial
## Min. :2021-04-27 Burrow :25 245-752: 5
## 1st Qu.:2021-04-28 Open :75 245-760: 5
## Median :2021-04-29 Shade - Full : 4 252-884: 5
## Mean :2021-04-29 Shade - Partial: 7 229-060: 4
## 3rd Qu.:2021-05-01 229-063: 4
## Max. :2021-05-02 229-069: 4
## (Other):84
## DateTime hour_of_day bpmTempC dwsTempOut
## Min. :2021-05-01 10:03:00.00 Min. : 2.00 Mode:logical Mode:logical
## 1st Qu.:2021-05-01 13:39:30.00 1st Qu.:10.50 NA's:111 NA's:111
## Median :2021-05-01 16:00:00.00 Median :13.00
## Mean :2021-05-01 21:58:14.11 Mean :12.08
## 3rd Qu.:2021-05-02 11:00:30.00 3rd Qu.:15.00
## Max. :2021-05-02 15:46:00.00 Max. :17.00
## NA's :60 NA's :60
## ModTempC
## Mode:logical
## NA's:111
##
##
##
##
## Now, put together the original and messy-now-clean dataframes, and add some formatting:
microhab <- microhab_clean %>%
rbind(microhab_early_all) %>%
dplyr::filter(complete.cases(microhabitat)) %>%
left_join(group_dates, by = 'date_only') %>%
left_join(liz_IDs_only, by = c('Serial' = 'radio_collar_serial')) %>%
mutate(above_below = factor(case_when(microhabitat == "Burrow" ~ "Belowground",
microhabitat != "Burrow" ~ "Aboveground"),
levels = c("Aboveground", "Belowground")))
summary(microhab)## Serial DateTime date_only
## 245-744: 100 Min. :2021-05-01 10:03:00.0 Min. :2021-04-27
## 245-741: 96 1st Qu.:2021-05-31 13:49:00.0 1st Qu.:2021-05-29
## 229-078: 93 Median :2021-06-22 13:29:00.0 Median :2021-06-21
## 229-079: 91 Mean :2021-06-15 06:38:36.9 Mean :2021-06-13
## 245-745: 91 3rd Qu.:2021-06-30 14:27:00.0 3rd Qu.:2021-06-30
## 245-748: 91 Max. :2021-07-13 13:27:00.0 Max. :2021-07-13
## (Other):1207 NA's :60
## hour_of_day microhabitat bpmTempC dwsTempOut
## Min. : 2 Open :540 Min. : 8.028 Min. :10.11
## 1st Qu.:10 Shade - Partial:147 1st Qu.:32.828 1st Qu.:26.61
## Median :12 Shade - Full :214 Median :36.267 Median :31.28
## Mean :12 Burrow :868 Mean :35.633 Mean :30.25
## 3rd Qu.:14 3rd Qu.:39.792 3rd Qu.:34.60
## Max. :18 Max. :49.823 Max. :40.58
## NA's :60 NA's :143 NA's :166
## ModTempC date_chunks individual_ID sex_M_F
## Min. :21.24 Apr 26 - May 6:159 M-09 : 100 Female: 535
## 1st Qu.:39.50 May 10 - 20 :158 F-10 : 96 Male :1234
## Median :44.48 May 21 - 31 :181 M-14 : 93
## Mean :43.77 Jun 1 - 11 :181 M-08 : 91
## 3rd Qu.:48.99 Jun 12 - 22 :253 M-11 : 91
## Max. :57.86 Jun 23 - Jul 3:579 M-19 : 91
## NA's :111 Jul 4 - 13 :258 (Other):1207
## above_below
## Aboveground:901
## Belowground:868
##
##
##
##
## Thermal Metrics
For each lizard, we want to calculate some assessment of how hot they let themselves be, and their thermoregulatory accuracy.
Daily Maximums
daily_max <- temp_clean %>%
# quick look to use
#dplyr::filter(date_chunks %in% c("Apr 26 - May 6", "May 10 - 20")) %>%
# first do daily maximums
group_by(date_only, date_chunks, Serial) %>%
mutate(Tb_percentile_50 = quantile(ppTempC, 0.5),
Tb_percentile_60 = quantile(ppTempC, 0.6),
Tb_percentile_70 = quantile(ppTempC, 0.7),
Tb_percentile_80 = quantile(ppTempC, 0.8),
Tb_percentile_90 = quantile(ppTempC, 0.9)) %>%
ungroup() %>%
# select only the data between the 80-90th percentiles of Tb for each lizard
dplyr::filter(ppTempC > Tb_percentile_80 & ppTempC < Tb_percentile_90) %>%
# calculate the mean Tb for that subset for each lizard per day
group_by(date_only, date_chunks, Serial) %>%
summarise(daily_Tb_percentile_80_90 = mean(ppTempC)) %>%
# then do the mean Tb between 80-90%ile for each lizard per time chunk (~11 days)
group_by(date_chunks, Serial) %>%
summarise(chunk_avg_Tb_percentile_80_90 = mean(daily_Tb_percentile_80_90))## `summarise()` has grouped output by 'date_only', 'date_chunks'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'date_chunks'. You can override using the
## `.groups` argument.summary(daily_max)## date_chunks Serial chunk_avg_Tb_percentile_80_90
## Apr 26 - May 6:35 229-065: 7 Min. :18.89
## May 10 - 20 :32 229-069: 7 1st Qu.:37.51
## May 21 - 31 :29 229-070: 7 Median :39.38
## Jun 1 - 11 :23 229-072: 7 Mean :38.37
## Jun 12 - 22 :19 229-078: 7 3rd Qu.:40.76
## Jun 23 - Jul 3:16 245-742: 7 Max. :46.85
## Jul 4 - 13 :13 (Other):125Db
From Kat Ivey’s 2020 paper, Tset range is: 32.3±1.2 - 37.5±1.1 deg C.
From Nicole Gaudenti’s 2021 paper, Tset range: 33.2-37.9 deg C.
I will use the widest margin created by both: 32.3 (Ivey) - 37.9 (Gaudenti)
Sunrise during our study: 6:05-6:40 am –> middle ~~6:20? Sunset during our study: 20:24-21:07 (8:24-9:07 pm) –> middle ~~20:40? (entire study during DST) from NOAA solar calculator
Or, just use the same daily time range as Ivey did? (0700-1900)
Although Db is usually taken as the absolute value, I will use pos/neg values, because I want to know the directionality of Db.
Db <- temp_clean %>%
# we only want daytime Tb for this
# currently using Kat's limits
dplyr::filter(hour_of_day <=19 & hour_of_day >=7) %>%
mutate(Tpref_low = 32.3,
Tpref_high = 37.9) %>%
# calculate Db for every Tb measurement
mutate(Db = case_when(ppTempC < Tpref_low ~ ppTempC - Tpref_low, # for pos/neg version
#ppTempC < Tpref_low ~ Tpref_low - ppTempC, # for abs value version
ppTempC > Tpref_high ~ ppTempC - Tpref_high,
ppTempC < Tpref_high & ppTempC > Tpref_low ~ 0)) %>%
# then avg Db by lizard & date chunk
group_by(date_chunks, Serial) %>%
summarise(chunk_mean_Db_deg_C_diff = mean(Db))## `summarise()` has grouped output by 'date_chunks'. You can override using the
## `.groups` argument.summary(Db)## date_chunks Serial chunk_mean_Db_deg_C_diff
## Apr 26 - May 6:35 229-065: 7 Min. :-16.4294
## May 10 - 20 :32 229-069: 7 1st Qu.: -0.2364
## May 21 - 31 :30 229-070: 7 Median : 0.6259
## Jun 1 - 11 :25 229-072: 7 Mean : -0.1235
## Jun 12 - 22 :21 229-078: 7 3rd Qu.: 1.3575
## Jun 23 - Jul 3:19 229-079: 7 Max. : 8.8909
## Jul 4 - 13 :18 (Other):138Microhabitat Use
Proportions of Observations
To calculate microhabitat use, I literally just calculate the proportion of observations for a given individual/time range that were in each microhabitat.
So, for each time chunk, for each lizard (radio-collar serial), I need to get the following values:
- number of total observations
- number of observations in each microhabitat
- proportion of observations in each microhabitat
- aboveground temps based on Ian’s model ?
We do not have a lot of obs for each lizard for every time chunk, so how many observations should there be to make the proportions valid/reasonable? One potential way to do this could be to at least have 1-2 obs for every microhabitat type we are trying to observe (so, minimum points per date chunk either 4 or 8, once I get rid of the irrelevant ones).
# where to cut off n obs?
microhab %>%
# get only 1 per lizard
group_by(date_chunks, Serial) %>%
summarise(n_obs = n()) %>%
ungroup() %>%
# How many lizards have each (n) obs?
group_by(n_obs) %>%
summarise(n = n())## `summarise()` has grouped output by 'date_chunks'. You can override using the
## `.groups` argument.## # A tibble: 26 × 2
## n_obs n
## <int> <int>
## 1 1 10
## 2 2 11
## 3 3 6
## 4 4 15
## 5 5 28
## 6 6 19
## 7 7 12
## 8 8 4
## 9 9 8
## 10 10 8
## # … with 16 more rowsIt seems like a good cut-off is a minimum of 4 observations per lizard, per time chunk. Otherwise, we would lose a LOT of data.
# micro use BY individual lizard
MH_use_by_indiv <- microhab %>%
group_by(Serial, date_chunks, microhabitat) %>%
summarise(n_obs_micro = n()) %>%
group_by(Serial, date_chunks) %>%
mutate(n_obs_total = sum(n_obs_micro)) %>%
ungroup() %>%
# expect that this takes out 42 obs/rows
dplyr::filter(n_obs_total >= 4) %>%
mutate(MH_use_proportion = n_obs_micro/n_obs_total) #%>%## `summarise()` has grouped output by 'Serial', 'date_chunks'. You can override
## using the `.groups` argument.# check that all proportions for a given lizard and time chunk = 1
# group_by(Serial, date_chunks) %>%
# summarise(prop_total = sum(MH_use_proportion))
head(MH_use_by_indiv)## # A tibble: 6 × 6
## Serial date_chunks microhabitat n_obs_micro n_obs_total MH_use_propor…¹
## <fct> <fct> <fct> <int> <int> <dbl>
## 1 229-059 Apr 26 - May 6 Open 3 4 0.75
## 2 229-059 Apr 26 - May 6 Burrow 1 4 0.25
## 3 229-060 Apr 26 - May 6 Open 4 6 0.667
## 4 229-060 Apr 26 - May 6 Burrow 2 6 0.333
## 5 229-060 May 10 - 20 Open 2 5 0.4
## 6 229-060 May 10 - 20 Shade - Partial 1 5 0.2
## # … with abbreviated variable name ¹MH_use_proportionsummary(MH_use_by_indiv)## Serial date_chunks microhabitat n_obs_micro
## 245-741 : 23 Apr 26 - May 6:60 Open :120 Min. : 1.00
## 229-079 : 22 May 10 - 20 :64 Shade - Partial: 74 1st Qu.: 1.00
## 245-745 : 22 May 21 - 31 :70 Shade - Full : 67 Median : 3.00
## 245-744 : 21 Jun 1 - 11 :59 Burrow :116 Mean : 4.56
## 245-762a: 19 Jun 12 - 22 :50 3rd Qu.: 5.00
## 245-761 : 18 Jun 23 - Jul 3:45 Max. :32.00
## (Other) :252 Jul 4 - 13 :29
## n_obs_total MH_use_proportion
## Min. : 4.00 Min. :0.0303
## 1st Qu.: 5.00 1st Qu.:0.1892
## Median : 9.00 Median :0.3077
## Mean :11.65 Mean :0.4005
## 3rd Qu.:14.00 3rd Qu.:0.6000
## Max. :37.00 Max. :1.0000
## # micro use for all lizards in the study, aggregated
MH_use_across_indivs <- microhab %>%
group_by(date_chunks, microhabitat) %>%
summarise(n_obs_micro = n(),
mean_mod_temp = mean(ModTempC, na.rm = T)) %>%
group_by(date_chunks) %>%
mutate(n_obs_total = sum(n_obs_micro),
mean_mod_temp2 = mean(mean_mod_temp, na.rm = T)) %>%
ungroup() %>%
mutate(MH_use_proportion = n_obs_micro/n_obs_total)## `summarise()` has grouped output by 'date_chunks'. You can override using the
## `.groups` argument.head(MH_use_across_indivs)## # A tibble: 6 × 7
## date_chunks microhabitat n_obs_micro mean_mod_…¹ n_obs…² mean_…³ MH_us…⁴
## <fct> <fct> <int> <dbl> <int> <dbl> <dbl>
## 1 Apr 26 - May 6 Open 106 37.3 159 40.5 0.667
## 2 Apr 26 - May 6 Shade - Partial 13 42.7 159 40.5 0.0818
## 3 Apr 26 - May 6 Shade - Full 8 41.4 159 40.5 0.0503
## 4 Apr 26 - May 6 Burrow 32 40.5 159 40.5 0.201
## 5 May 10 - 20 Open 82 42.3 158 44.3 0.519
## 6 May 10 - 20 Shade - Partial 23 44.3 158 44.3 0.146
## # … with abbreviated variable names ¹mean_mod_temp, ²n_obs_total,
## # ³mean_mod_temp2, ⁴MH_use_proportion# micro use for all lizards, BY SEX
MH_use_across_indivs_sex <- microhab %>%
group_by(date_chunks, microhabitat, sex_M_F) %>%
summarise(n_obs_micro = n(),
mean_mod_temp = mean(ModTempC, na.rm = T)) %>%
group_by(date_chunks, sex_M_F) %>%
mutate(n_obs_total = sum(n_obs_micro),
mean_mod_temp2 = mean(mean_mod_temp, na.rm = T)) %>%
ungroup() %>%
mutate(MH_use_proportion = n_obs_micro/n_obs_total)## `summarise()` has grouped output by 'date_chunks', 'microhabitat'. You can
## override using the `.groups` argument.summary(MH_use_across_indivs_sex)## date_chunks microhabitat sex_M_F n_obs_micro
## Apr 26 - May 6:8 Open :14 Female:28 Min. : 1.00
## May 10 - 20 :8 Shade - Partial:14 Male :28 1st Qu.: 9.00
## May 21 - 31 :8 Shade - Full :14 Median : 16.00
## Jun 1 - 11 :8 Burrow :14 Mean : 31.59
## Jun 12 - 22 :8 3rd Qu.: 32.50
## Jun 23 - Jul 3:8 Max. :269.00
## Jul 4 - 13 :8
## mean_mod_temp n_obs_total mean_mod_temp2 MH_use_proportion
## Min. :33.22 Min. : 52.0 Min. :39.76 Min. :0.01724
## 1st Qu.:40.21 1st Qu.: 70.0 1st Qu.:40.56 1st Qu.:0.09049
## Median :43.50 Median : 85.0 Median :43.75 Median :0.15281
## Mean :43.08 Mean :126.4 Mean :43.08 Mean :0.25000
## 3rd Qu.:45.85 3rd Qu.:141.0 3rd Qu.:44.92 3rd Qu.:0.37761
## Max. :49.76 Max. :438.0 Max. :45.21 Max. :0.77000
## calculate the proportion of time aboveground:
prop_above_by_indiv <- MH_use_by_indiv %>%
dplyr::filter(microhabitat == "Burrow") %>%
mutate(prop_above = 1 - MH_use_proportion)
prop_above_across_indiv <- MH_use_across_indivs %>%
dplyr::filter(microhabitat == "Burrow") %>%
mutate(prop_above = 1 - MH_use_proportion)
prop_above_across_indiv_by_sex <- MH_use_across_indivs_sex %>%
dplyr::filter(microhabitat == "Burrow") %>%
mutate(prop_above = 1 - MH_use_proportion)Onset of Estivation
estivating <- microhab %>%
# for each lizard
group_by(Serial) %>%
# what was the latest date we tracked them?
mutate(date_last_obs = max(date_only)) %>%
ungroup() %>%
# now, subset aboveground obs only
dplyr::filter(microhabitat != "Burrow") %>%
# what was the latest date we observed each lizard ABOVEground?
group_by(Serial) %>%
mutate(date_last_aboveground_obs = max(date_only),
est_obs = date_last_obs - date_last_aboveground_obs) %>%
ungroup() %>%
# make sure we tracked them after "last" aboveground observation
dplyr::filter(date_last_aboveground_obs < date_last_obs) %>%
# for each lizard
group_by(Serial) %>%
summarise(est_onset = mean(date_last_aboveground_obs),
est_obs_days = as.numeric(mean(est_obs))
) %>%
dplyr::filter(est_obs_days >1)
summary(estivating)## Serial est_onset est_obs_days
## 229-059 : 1 Min. :2021-05-01 Min. : 3.00
## 229-060 : 1 1st Qu.:2021-05-20 1st Qu.: 7.00
## 229-063 : 1 Median :2021-06-11 Median :18.00
## 229-065 : 1 Mean :2021-06-07 Mean :19.05
## 229-069 : 1 3rd Qu.:2021-06-25 3rd Qu.:27.00
## 229-070a: 1 Max. :2021-07-08 Max. :63.00
## (Other) :15Max Temps ~ Hydric
Prep Data
max_temps_hydric <- liz_dat_apr_may_only %>%
left_join(daily_max, by = c('date_chunks',
'radio_collar_serial' = 'Serial')) %>%
# remove NAs
dplyr::filter(complete.cases(chunk_avg_Tb_percentile_80_90)) %>%
# remove the one outlier wayyyyy below the other values
dplyr::filter(chunk_avg_Tb_percentile_80_90>30)
summary(max_temps_hydric)## capture_date_time capture_date radio_collar_serial
## Min. :2021-04-23 10:39:00 Min. :2021-04-23 229-059: 2
## 1st Qu.:2021-04-23 13:51:15 1st Qu.:2021-04-23 229-060: 2
## Median :2021-04-24 10:45:00 Median :2021-04-24 229-063: 2
## Mean :2021-04-28 01:01:05 Mean :2021-04-28 229-069: 2
## 3rd Qu.:2021-05-07 11:03:45 3rd Qu.:2021-05-07 229-072: 2
## Max. :2021-05-08 12:27:00 Max. :2021-05-08 229-077: 2
## NA's :4 (Other):40
## PIT_tag_ID individual_ID sex_M_F mass_g SVL_mm
## Length:52 F-01 : 2 Female:26 Min. :26.00 Min. : 85.00
## Class :character F-05 : 2 Male :26 1st Qu.:33.00 1st Qu.: 94.75
## Mode :character F-06 : 2 Median :37.00 Median : 99.00
## F-11 : 2 Mean :38.72 Mean :100.19
## F-12 : 2 3rd Qu.:43.45 3rd Qu.:107.25
## F-13 : 2 Max. :53.60 Max. :122.00
## (Other):40
## tmt tmt_mass_change_g capture_to_msmt hematocrit_percent
## Water Tmt:23 Min. :-3.0000 Min. : 11.00 Min. :23.00
## Sham Tmt :29 1st Qu.:-1.0000 1st Qu.: 61.25 1st Qu.:30.00
## No Tmt : 0 Median : 0.0000 Median : 152.00 Median :36.00
## Mean :-0.3423 Mean : 398.88 Mean :34.51
## 3rd Qu.: 0.0000 3rd Qu.: 263.00 3rd Qu.:38.00
## Max. : 3.0000 Max. :1665.00 Max. :58.00
## NA's :4 NA's :1
## cloacal_temp_C CEWL_g_m2h msmt_temp_C msmt_RH_percent
## Min. :24.00 Min. : 0.650 Min. :19.07 Min. :12.20
## 1st Qu.:32.00 1st Qu.: 7.677 1st Qu.:28.68 1st Qu.:14.53
## Median :33.00 Median : 9.262 Median :30.62 Median :17.87
## Mean :32.44 Mean : 9.570 Mean :29.44 Mean :19.34
## 3rd Qu.:34.00 3rd Qu.:11.723 3rd Qu.:31.88 3rd Qu.:22.23
## Max. :36.00 Max. :21.673 Max. :33.55 Max. :35.83
##
## e_s_kPa msmt_VPD_kPa osmolality_mmol_kg_mean month week
## Min. :2.205 Min. :1.415 Min. :316.0 April:34 1 :34
## 1st Qu.:3.931 1st Qu.:3.056 1st Qu.:346.6 May :18 3 :18
## Median :4.396 Median :3.694 Median :363.5 July : 0 12: 0
## Mean :4.178 Mean :3.412 Mean :365.9
## 3rd Qu.:4.721 3rd Qu.:4.037 3rd Qu.:380.2
## Max. :5.188 Max. :4.319 Max. :436.5
## NA's :1
## week_num ultrasound_date gravid_Y_N number_eggs
## Min. :1.000 Min. :2021-04-24 Gravid :13 Min. :3.000
## 1st Qu.:1.000 1st Qu.:2021-04-24 Not Gravid:11 1st Qu.:3.000
## Median :1.000 Median :2021-04-24 NA's :28 Median :4.000
## Mean :1.692 Mean :2021-04-29 Mean :3.692
## 3rd Qu.:3.000 3rd Qu.:2021-05-08 3rd Qu.:4.000
## Max. :3.000 Max. :2021-05-08 Max. :4.000
## NA's :28 NA's :39
## largest_egg_diameter_cm largest_egg_circumference_cm stage
## Min. :0.640 Min. :1.820 small round: 1
## 1st Qu.:0.870 1st Qu.:2.450 large round: 7
## Median :0.990 Median :2.890 soft : 2
## Mean :1.118 Mean :3.171 soft oblong: 0
## 3rd Qu.:1.400 3rd Qu.:3.970 firm oblong: 2
## Max. :1.800 Max. :4.740 hard oblong: 1
## NA's :39 NA's :39 NA's :39
## dev_point SMI date_chunks
## Min. :1.000 Min. :27.04 Apr 26 - May 6:34
## 1st Qu.:2.000 1st Qu.:33.00 May 10 - 20 :18
## Median :2.000 Median :34.61 May 21 - 31 : 0
## Mean :2.769 Mean :35.48 Jun 1 - 11 : 0
## 3rd Qu.:3.000 3rd Qu.:38.49 Jun 12 - 22 : 0
## Max. :5.000 Max. :45.63 Jun 23 - Jul 3: 0
## NA's :39 Jul 4 - 13 : 0
## chunk_avg_Tb_percentile_80_90
## Min. :35.53
## 1st Qu.:38.36
## Median :39.14
## Mean :39.06
## 3rd Qu.:39.73
## Max. :41.64
## Plots
ggplot(max_temps_hydric) +
geom_point(aes(x = CEWL_g_m2h,
y = chunk_avg_Tb_percentile_80_90,
color = date_chunks)) +
geom_smooth(aes(x = CEWL_g_m2h,
y = chunk_avg_Tb_percentile_80_90),
se = F,
method = "lm",
formula = y~x)
ggplot(max_temps_hydric) +
geom_point(aes(x = osmolality_mmol_kg_mean,
y = chunk_avg_Tb_percentile_80_90,
color = date_chunks)) +
geom_smooth(aes(x = osmolality_mmol_kg_mean,
y = chunk_avg_Tb_percentile_80_90),
se = F,
method = "lm",
formula = y~x)## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 1 rows containing missing values (`geom_point()`).
Models
# CEWL
max_temp_CEWL_mod <- lmerTest::lmer(data = max_temps_hydric,
chunk_avg_Tb_percentile_80_90 ~ CEWL_g_m2h +
(1|individual_ID))
summary(max_temp_CEWL_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: chunk_avg_Tb_percentile_80_90 ~ CEWL_g_m2h + (1 | individual_ID)
## Data: max_temps_hydric
##
## REML criterion at convergence: 167.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.06899 -0.46408 -0.02622 0.47948 1.67363
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 0.6833 0.8266
## Residual 0.7950 0.8916
## Number of obs: 52, groups: individual_ID, 34
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 39.17879 0.39372 47.13737 99.509 <2e-16 ***
## CEWL_g_m2h -0.02380 0.03627 37.66944 -0.656 0.516
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## CEWL_g_m2h -0.875anova(max_temp_CEWL_mod)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## CEWL_g_m2h 0.34231 0.34231 1 37.669 0.4306 0.5157plot(max_temp_CEWL_mod)
# osmolality
max_temp_osml_mod <- lmerTest::lmer(data = max_temps_hydric,
chunk_avg_Tb_percentile_80_90 ~ osmolality_mmol_kg_mean +
(1|individual_ID))
summary(max_temp_osml_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: chunk_avg_Tb_percentile_80_90 ~ osmolality_mmol_kg_mean + (1 |
## individual_ID)
## Data: max_temps_hydric
##
## REML criterion at convergence: 165.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.09186 -0.40928 0.00062 0.49188 1.73478
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 0.7427 0.8618
## Residual 0.7014 0.8375
## Number of obs: 51, groups: individual_ID, 33
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 35.097457 2.124721 43.207436 16.519 <2e-16 ***
## osmolality_mmol_kg_mean 0.010601 0.005795 43.112375 1.829 0.0743 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## osmllty_m__ -0.996anova(max_temp_osml_mod)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## osmolality_mmol_kg_mean 2.3473 2.3473 1 43.112 3.3468 0.07426 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(max_temp_osml_mod)
# mass
max_temp_mass_mod <- lmerTest::lmer(data = max_temps_hydric,
chunk_avg_Tb_percentile_80_90 ~ mass_g +
(1|individual_ID))
summary(max_temp_mass_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: chunk_avg_Tb_percentile_80_90 ~ mass_g + (1 | individual_ID)
## Data: max_temps_hydric
##
## REML criterion at convergence: 167.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.18341 -0.45928 0.00822 0.44081 1.80182
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 0.6609 0.8130
## Residual 0.7805 0.8834
## Number of obs: 52, groups: individual_ID, 34
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 37.80830 0.91545 38.36858 41.300 <2e-16 ***
## mass_g 0.02985 0.02335 37.68275 1.278 0.209
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## mass_g -0.979anova(max_temp_mass_mod)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## mass_g 1.2754 1.2754 1 37.683 1.6342 0.2089plot(max_temp_mass_mod)
Db ~ Hydric
Prep Data
db_hydric <- liz_dat_apr_may_only %>%
left_join(Db, by = c('date_chunks',
'radio_collar_serial' = 'Serial')) %>%
# remove NAs
dplyr::filter(complete.cases(chunk_mean_Db_deg_C_diff)) %>%
# remove outliers
dplyr::filter(chunk_mean_Db_deg_C_diff > -4)
summary(db_hydric)## capture_date_time capture_date radio_collar_serial
## Min. :2021-04-23 10:39:00.00 Min. :2021-04-23 229-059: 2
## 1st Qu.:2021-04-23 13:42:30.00 1st Qu.:2021-04-23 229-060: 2
## Median :2021-04-24 10:49:00.00 Median :2021-04-24 229-063: 2
## Mean :2021-04-28 03:16:33.19 Mean :2021-04-28 229-069: 2
## 3rd Qu.:2021-05-07 11:07:30.00 3rd Qu.:2021-05-07 229-072: 2
## Max. :2021-05-08 12:27:00.00 Max. :2021-05-08 229-077: 2
## NA's :4 (Other):39
## PIT_tag_ID individual_ID sex_M_F mass_g SVL_mm
## Length:51 F-01 : 2 Female:26 Min. :26.00 Min. : 85.0
## Class :character F-05 : 2 Male :25 1st Qu.:33.00 1st Qu.: 94.5
## Mode :character F-06 : 2 Median :37.00 Median : 99.0
## F-11 : 2 Mean :38.79 Mean :100.2
## F-12 : 2 3rd Qu.:43.60 3rd Qu.:107.5
## F-13 : 2 Max. :53.60 Max. :122.0
## (Other):39
## tmt tmt_mass_change_g capture_to_msmt hematocrit_percent
## Water Tmt:22 Min. :-3.000 Min. : 11.0 Min. :23.00
## Sham Tmt :29 1st Qu.:-1.000 1st Qu.: 60.5 1st Qu.:30.00
## No Tmt : 0 Median : 0.000 Median : 150.0 Median :36.00
## Mean :-0.349 Mean : 372.0 Mean :34.64
## 3rd Qu.: 0.000 3rd Qu.: 259.0 3rd Qu.:38.00
## Max. : 3.000 Max. :1665.0 Max. :58.00
## NA's :4 NA's :1
## cloacal_temp_C CEWL_g_m2h msmt_temp_C msmt_RH_percent
## Min. :24.00 Min. : 0.650 Min. :19.07 Min. :12.20
## 1st Qu.:32.00 1st Qu.: 7.377 1st Qu.:28.74 1st Qu.:14.53
## Median :33.00 Median : 9.207 Median :30.68 Median :17.07
## Mean :32.51 Mean : 9.574 Mean :29.61 Mean :19.08
## 3rd Qu.:34.00 3rd Qu.:11.886 3rd Qu.:31.88 3rd Qu.:22.12
## Max. :36.00 Max. :21.673 Max. :33.55 Max. :35.83
##
## e_s_kPa msmt_VPD_kPa osmolality_mmol_kg_mean month week
## Min. :2.205 Min. :1.415 Min. :316.0 April:33 1 :33
## 1st Qu.:3.945 1st Qu.:3.068 1st Qu.:346.0 May :18 3 :18
## Median :4.411 Median :3.703 Median :362.8 July : 0 12: 0
## Mean :4.211 Mean :3.446 Mean :365.3
## 3rd Qu.:4.723 3rd Qu.:4.038 3rd Qu.:379.2
## Max. :5.188 Max. :4.319 Max. :436.5
## NA's :1
## week_num ultrasound_date gravid_Y_N number_eggs
## Min. :1.000 Min. :2021-04-24 Gravid :13 Min. :3.000
## 1st Qu.:1.000 1st Qu.:2021-04-24 Not Gravid:11 1st Qu.:3.000
## Median :1.000 Median :2021-04-24 NA's :27 Median :4.000
## Mean :1.706 Mean :2021-04-29 Mean :3.692
## 3rd Qu.:3.000 3rd Qu.:2021-05-08 3rd Qu.:4.000
## Max. :3.000 Max. :2021-05-08 Max. :4.000
## NA's :27 NA's :38
## largest_egg_diameter_cm largest_egg_circumference_cm stage
## Min. :0.640 Min. :1.820 small round: 1
## 1st Qu.:0.870 1st Qu.:2.450 large round: 7
## Median :0.990 Median :2.890 soft : 2
## Mean :1.118 Mean :3.171 soft oblong: 0
## 3rd Qu.:1.400 3rd Qu.:3.970 firm oblong: 2
## Max. :1.800 Max. :4.740 hard oblong: 1
## NA's :38 NA's :38 NA's :38
## dev_point SMI date_chunks chunk_mean_Db_deg_C_diff
## Min. :1.000 Min. :27.04 Apr 26 - May 6:33 Min. :-0.5788
## 1st Qu.:2.000 1st Qu.:33.00 May 10 - 20 :18 1st Qu.: 0.1958
## Median :2.000 Median :34.71 May 21 - 31 : 0 Median : 0.5853
## Mean :2.769 Mean :35.51 Jun 1 - 11 : 0 Mean : 0.6019
## 3rd Qu.:3.000 3rd Qu.:38.52 Jun 12 - 22 : 0 3rd Qu.: 0.8768
## Max. :5.000 Max. :45.63 Jun 23 - Jul 3: 0 Max. : 1.7634
## NA's :38 Jul 4 - 13 : 0Plots
ggplot(db_hydric) +
geom_point(aes(x = CEWL_g_m2h,
y = chunk_mean_Db_deg_C_diff,
color = date_chunks)) +
geom_smooth(aes(x = CEWL_g_m2h,
y = chunk_mean_Db_deg_C_diff),
se = F,
method = "lm",
formula = y~x)
ggplot(db_hydric) +
geom_point(aes(x = osmolality_mmol_kg_mean,
y = chunk_mean_Db_deg_C_diff,
color = date_chunks)) +
geom_smooth(aes(x = osmolality_mmol_kg_mean,
y = chunk_mean_Db_deg_C_diff),
se = F,
method = "lm",
formula = y~x)## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 1 rows containing missing values (`geom_point()`).
Models
# CEWL
db_CEWL_mod <- lmerTest::lmer(data = db_hydric,
chunk_mean_Db_deg_C_diff ~ CEWL_g_m2h +
(1|individual_ID))
summary(db_CEWL_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: chunk_mean_Db_deg_C_diff ~ CEWL_g_m2h + (1 | individual_ID)
## Data: db_hydric
##
## REML criterion at convergence: 93.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.62121 -0.65120 -0.01753 0.38077 1.92640
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 0.08578 0.2929
## Residual 0.23890 0.4888
## Number of obs: 51, groups: individual_ID, 33
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.73930 0.19320 48.78761 3.827 0.00037 ***
## CEWL_g_m2h -0.01774 0.01813 45.65577 -0.978 0.33304
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## CEWL_g_m2h -0.894anova(db_CEWL_mod)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## CEWL_g_m2h 0.22869 0.22869 1 45.656 0.9572 0.333plot(db_CEWL_mod)
# osmolality
db_osml_mod <- lmerTest::lmer(data = db_hydric,
chunk_mean_Db_deg_C_diff ~ osmolality_mmol_kg_mean +
(1|individual_ID))
summary(db_osml_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## chunk_mean_Db_deg_C_diff ~ osmolality_mmol_kg_mean + (1 | individual_ID)
## Data: db_hydric
##
## REML criterion at convergence: 92.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.63855 -0.61158 -0.09174 0.58990 1.76793
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 0.09429 0.3071
## Residual 0.21515 0.4638
## Number of obs: 50, groups: individual_ID, 32
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.391181 1.054351 47.767585 -1.319 0.1933
## osmolality_mmol_kg_mean 0.005411 0.002881 47.802411 1.878 0.0665 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## osmllty_m__ -0.997anova(db_osml_mod)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## osmolality_mmol_kg_mean 0.7589 0.7589 1 47.802 3.5274 0.06647 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(db_osml_mod)
# mass
db_mass_mod <- lmerTest::lmer(data = db_hydric,
chunk_mean_Db_deg_C_diff ~ mass_g +
(1|individual_ID))
summary(db_mass_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: chunk_mean_Db_deg_C_diff ~ mass_g + (1 | individual_ID)
## Data: db_hydric
##
## REML criterion at convergence: 93.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.76058 -0.64930 -0.09152 0.48856 1.90690
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 0.07312 0.2704
## Residual 0.24566 0.4956
## Number of obs: 51, groups: individual_ID, 33
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.08040 0.42216 35.83347 0.190 0.850
## mass_g 0.01279 0.01071 34.78071 1.194 0.241
##
## Correlation of Fixed Effects:
## (Intr)
## mass_g -0.980anova(db_mass_mod)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## mass_g 0.35018 0.35018 1 34.781 1.4255 0.2406plot(db_mass_mod)
Microhabitat ~ Hydric
Prep Data
I actually need to prep the microhabitat data- the metric I will investigate is proportion of time aboveground.
prop_above_hydric <- prop_above_by_indiv %>%
dplyr::filter(date_chunks %in% c("Apr 26 - May 6", "May 10 - 20"))micro_hydric <- liz_dat_apr_may_only %>%
left_join(prop_above_hydric, by = c('date_chunks',
'radio_collar_serial' = 'Serial')) %>%
# remove NAs
dplyr::filter(complete.cases(prop_above))
summary(micro_hydric)## capture_date_time capture_date radio_collar_serial
## Min. :2021-04-23 10:43:00.00 Min. :2021-04-23 229-060: 2
## 1st Qu.:2021-04-23 13:09:00.00 1st Qu.:2021-04-23 252-884: 2
## Median :2021-04-24 10:23:00.00 Median :2021-04-24 252-887: 2
## Mean :2021-04-27 19:59:08.18 Mean :2021-04-28 229-059: 1
## 3rd Qu.:2021-05-04 12:35:45.00 3rd Qu.:2021-05-07 229-066: 1
## Max. :2021-05-08 12:27:00.00 Max. :2021-05-08 229-070: 1
## NA's :2 (Other):15
## PIT_tag_ID individual_ID sex_M_F mass_g SVL_mm
## Length:24 F-11 : 2 Female:14 Min. :26.00 Min. : 85.00
## Class :character F-17 : 2 Male :10 1st Qu.:32.00 1st Qu.: 94.00
## Mode :character M-04 : 2 Median :35.30 Median : 96.50
## F-01 : 1 Mean :36.76 Mean : 98.46
## F-03 : 1 3rd Qu.:38.45 3rd Qu.:104.25
## F-05 : 1 Max. :52.70 Max. :114.00
## (Other):15
## tmt tmt_mass_change_g capture_to_msmt hematocrit_percent
## Water Tmt:10 Min. :-3.0000 Min. : 25.0 Min. :23.00
## Sham Tmt :14 1st Qu.:-1.0000 1st Qu.: 77.5 1st Qu.:27.50
## No Tmt : 0 Median :-0.1500 Median : 157.0 Median :32.00
## Mean :-0.4833 Mean : 455.3 Mean :34.17
## 3rd Qu.: 0.0000 3rd Qu.: 270.2 3rd Qu.:40.00
## Max. : 1.8000 Max. :1662.0 Max. :58.00
## NA's :2 NA's :1
## cloacal_temp_C CEWL_g_m2h msmt_temp_C msmt_RH_percent
## Min. :24.00 Min. : 0.650 Min. :19.07 Min. :13.12
## 1st Qu.:32.00 1st Qu.: 8.443 1st Qu.:28.39 1st Qu.:14.66
## Median :33.00 Median : 9.430 Median :30.62 Median :16.68
## Mean :32.58 Mean :10.401 Mean :29.01 Mean :19.71
## 3rd Qu.:34.00 3rd Qu.:12.515 3rd Qu.:31.82 3rd Qu.:22.39
## Max. :36.00 Max. :21.673 Max. :33.07 Max. :35.83
##
## e_s_kPa msmt_VPD_kPa osmolality_mmol_kg_mean month week
## Min. :2.205 Min. :1.415 Min. :316.0 April:16 1 :16
## 1st Qu.:3.866 1st Qu.:2.987 1st Qu.:339.9 May : 8 3 : 8
## Median :4.396 Median :3.694 Median :354.5 July : 0 12: 0
## Mean :4.090 Mean :3.334 Mean :357.9
## 3rd Qu.:4.705 3rd Qu.:4.012 3rd Qu.:376.9
## Max. :5.049 Max. :4.309 Max. :402.0
## NA's :1
## week_num ultrasound_date gravid_Y_N number_eggs
## Min. :1.000 Min. :2021-04-24 Gravid : 8 Min. :3.000
## 1st Qu.:1.000 1st Qu.:2021-04-24 Not Gravid: 5 1st Qu.:3.000
## Median :1.000 Median :2021-04-24 NA's :11 Median :4.000
## Mean :1.667 Mean :2021-04-30 Mean :3.625
## 3rd Qu.:3.000 3rd Qu.:2021-05-08 3rd Qu.:4.000
## Max. :3.000 Max. :2021-05-08 Max. :4.000
## NA's :11 NA's :16
## largest_egg_diameter_cm largest_egg_circumference_cm stage
## Min. :0.6400 Min. :1.820 small round: 1
## 1st Qu.:0.8425 1st Qu.:2.385 large round: 5
## Median :0.9650 Median :2.885 soft : 1
## Mean :1.0488 Mean :2.976 soft oblong: 0
## 3rd Qu.:1.2700 3rd Qu.:3.618 firm oblong: 1
## Max. :1.5100 Max. :4.070 hard oblong: 0
## NA's :16 NA's :16 NA's :16
## dev_point SMI date_chunks microhabitat
## Min. :1.000 Min. :28.35 Apr 26 - May 6:16 Open : 0
## 1st Qu.:2.000 1st Qu.:32.48 May 10 - 20 : 8 Shade - Partial: 0
## Median :2.000 Median :34.28 May 21 - 31 : 0 Shade - Full : 0
## Mean :2.375 Mean :35.05 Jun 1 - 11 : 0 Burrow :24
## 3rd Qu.:2.250 3rd Qu.:35.75 Jun 12 - 22 : 0
## Max. :5.000 Max. :45.63 Jun 23 - Jul 3: 0
## NA's :16 Jul 4 - 13 : 0
## n_obs_micro n_obs_total MH_use_proportion prop_above
## Min. :1.000 Min. :4.000 Min. :0.1429 Min. :0.2500
## 1st Qu.:1.000 1st Qu.:4.750 1st Qu.:0.2000 1st Qu.:0.6000
## Median :1.000 Median :5.000 Median :0.2500 Median :0.7500
## Mean :1.708 Mean :5.292 Mean :0.3278 Mean :0.6722
## 3rd Qu.:2.000 3rd Qu.:6.000 3rd Qu.:0.4000 3rd Qu.:0.8000
## Max. :4.000 Max. :7.000 Max. :0.7500 Max. :0.8571
## Plots
ggplot(MH_use_by_indiv) +
geom_bar(aes(x = date_chunks,
y = MH_use_proportion,
fill = microhabitat
),
position = "fill",
stat = "identity") +
facet_wrap(~Serial) +
theme_classic()
ggplot(micro_hydric) +
geom_point(aes(x = CEWL_g_m2h,
y = prop_above,
color = date_chunks)) +
geom_smooth(aes(x = CEWL_g_m2h,
y = prop_above),
se = F,
method = "lm",
formula = y~x)
ggplot(micro_hydric) +
geom_point(aes(x = osmolality_mmol_kg_mean,
y = prop_above,
color = date_chunks)) +
geom_smooth(aes(x = osmolality_mmol_kg_mean,
y = prop_above),
se = F,
method = "lm",
formula = y~x)## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 1 rows containing missing values (`geom_point()`).
Models
# CEWL
micro_CEWL_mod <- lmerTest::lmer(data = micro_hydric,
prop_above ~ CEWL_g_m2h +
(1|individual_ID))
summary(micro_CEWL_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: prop_above ~ CEWL_g_m2h + (1 | individual_ID)
## Data: micro_hydric
##
## REML criterion at convergence: -6.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.9938 -0.2949 0.2102 0.3160 0.7489
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 0.026542 0.16292
## Residual 0.005813 0.07624
## Number of obs: 24, groups: individual_ID, 21
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.715491 0.084402 10.138401 8.477 6.46e-06 ***
## CEWL_g_m2h -0.004600 0.007328 7.028203 -0.628 0.55
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## CEWL_g_m2h -0.887anova(micro_CEWL_mod)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## CEWL_g_m2h 0.0022909 0.0022909 1 7.0282 0.3941 0.55plot(micro_CEWL_mod)
# osmolality
micro_osml_mod <- lmerTest::lmer(data = micro_hydric,
prop_above ~ osmolality_mmol_kg_mean +
(1|individual_ID))
summary(micro_osml_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: prop_above ~ osmolality_mmol_kg_mean + (1 | individual_ID)
## Data: micro_hydric
##
## REML criterion at convergence: -2.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.9601 -0.4517 0.1643 0.3356 0.8348
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 0.027168 0.16483
## Residual 0.005605 0.07487
## Number of obs: 23, groups: individual_ID, 20
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.9140781 0.4283049 7.1188832 2.134 0.0696 .
## osmolality_mmol_kg_mean -0.0007087 0.0011911 7.0127857 -0.595 0.5705
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## osmllty_m__ -0.996anova(micro_osml_mod)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## osmolality_mmol_kg_mean 0.0019848 0.0019848 1 7.0128 0.3541 0.5705plot(micro_osml_mod)
# mass
micro_mass_mod <- lmerTest::lmer(data = micro_hydric,
prop_above ~ mass_g +
(1|individual_ID))
summary(micro_mass_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: prop_above ~ mass_g + (1 | individual_ID)
## Data: micro_hydric
##
## REML criterion at convergence: -6.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.07862 -0.22959 0.07025 0.35562 1.04268
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 0.02671 0.16344
## Residual 0.00531 0.07287
## Number of obs: 24, groups: individual_ID, 21
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.502653 0.186882 20.318169 2.690 0.014 *
## mass_g 0.004510 0.004973 20.386310 0.907 0.375
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## mass_g -0.978anova(micro_mass_mod)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## mass_g 0.0043682 0.0043682 1 20.386 0.8226 0.375plot(micro_mass_mod)
Estivation ~ Hydric
Prep Data
est_hydric <- liz_dat_apr_may_only %>%
group_by(radio_collar_serial) %>%
summarise(mean_mass = mean(mass_g),
mean_SMI = mean(SMI),
mean_CEWL = mean(CEWL_g_m2h),
mean_osml = mean(osmolality_mmol_kg_mean, na.rm = T),
) %>%
left_join(estivating, by = c('radio_collar_serial' = 'Serial')) %>%
# remove NAs
dplyr::filter(complete.cases(est_onset)) %>%
mutate(est_onset_scaled = scale(as.numeric(est_onset)))
summary(est_hydric)## radio_collar_serial mean_mass mean_SMI mean_CEWL
## 229-059 : 1 Min. :28.00 Min. :27.22 Min. : 0.650
## 229-060 : 1 1st Qu.:36.95 1st Qu.:33.91 1st Qu.: 6.365
## 229-063 : 1 Median :38.90 Median :34.82 Median : 9.495
## 229-065 : 1 Mean :40.78 Mean :35.80 Mean : 9.339
## 229-069 : 1 3rd Qu.:46.00 3rd Qu.:37.68 3rd Qu.:11.352
## 229-070a: 1 Max. :55.40 Max. :42.58 Max. :18.743
## (Other) :15
## mean_osml est_onset est_obs_days est_onset_scaled.V1
## Min. :335.2 Min. :2021-05-01 Min. : 3.00 Min. :-1.7647589
## 1st Qu.:347.6 1st Qu.:2021-05-20 1st Qu.: 7.00 1st Qu.:-0.8723143
## Median :354.4 Median :2021-06-11 Median :18.00 Median : 0.1610426
## Mean :364.3 Mean :2021-06-07 Mean :19.05 Mean : 0.0000000
## 3rd Qu.:377.6 3rd Qu.:2021-06-25 3rd Qu.:27.00 3rd Qu.: 0.8186334
## Max. :418.5 Max. :2021-07-08 Max. :63.00 Max. : 1.4292534
## NA's :1also look at just July:
est_hydric_July <- liz_dat_july_only %>%
left_join(estivating, by = c('radio_collar_serial' = 'Serial')) %>%
# remove NAs
dplyr::filter(complete.cases(est_onset)) %>%
mutate(est_onset_scaled = scale(as.numeric(est_onset)))Plots
ggplot(est_hydric) +
geom_point(aes(x = mean_CEWL,
y = est_onset_scaled)) +
geom_smooth(aes(x = mean_CEWL,
y = est_onset_scaled),
se = F,
method = "lm",
formula = y~x)
ggplot(est_hydric) +
geom_point(aes(x = mean_osml,
y = est_onset_scaled)) +
geom_smooth(aes(x = mean_osml,
y = est_onset_scaled),
se = F,
method = "lm",
formula = y~x)## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 1 rows containing missing values (`geom_point()`).
Models
# CEWL
est_CEWL_mod <- lm(data = est_hydric, est_onset_scaled ~ mean_CEWL)
summary(est_CEWL_mod)##
## Call:
## lm(formula = est_onset_scaled ~ mean_CEWL, data = est_hydric)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.77722 -0.35069 0.01434 0.67241 1.49800
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.88047 0.51351 1.715 0.1027
## mean_CEWL -0.09428 0.05038 -1.871 0.0768 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9428 on 19 degrees of freedom
## Multiple R-squared: 0.1556, Adjusted R-squared: 0.1112
## F-statistic: 3.502 on 1 and 19 DF, p-value: 0.07677anova(est_CEWL_mod)## Analysis of Variance Table
##
## Response: est_onset_scaled
## Df Sum Sq Mean Sq F value Pr(>F)
## mean_CEWL 1 3.1126 3.11261 3.502 0.07677 .
## Residuals 19 16.8874 0.88881
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(est_CEWL_mod)
# osmolality
est_osml_mod <- lm(data = est_hydric, est_onset_scaled ~ mean_osml)
summary(est_osml_mod)##
## Call:
## lm(formula = est_onset_scaled ~ mean_osml, data = est_hydric)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.10255 -0.60198 0.09908 0.62680 1.22934
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.644281 3.187057 -2.712 0.0143 *
## mean_osml 0.023575 0.008732 2.700 0.0147 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8606 on 18 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.2883, Adjusted R-squared: 0.2487
## F-statistic: 7.29 on 1 and 18 DF, p-value: 0.01465anova(est_osml_mod)## Analysis of Variance Table
##
## Response: est_onset_scaled
## Df Sum Sq Mean Sq F value Pr(>F)
## mean_osml 1 5.3985 5.3985 7.2898 0.01465 *
## Residuals 18 13.3299 0.7406
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(est_osml_mod)
# osmolality JULY
est_osml_JULY_mod <- lm(data = est_hydric_July, est_onset_scaled ~ osmolality_mmol_kg_mean)
summary(est_osml_JULY_mod)##
## Call:
## lm(formula = est_onset_scaled ~ osmolality_mmol_kg_mean, data = est_hydric_July)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.22743 -0.79743 -0.03738 0.78900 1.61522
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.039091 7.080175 0.429 0.678
## osmolality_mmol_kg_mean -0.008547 0.019893 -0.430 0.678
##
## Residual standard error: 1.043 on 9 degrees of freedom
## Multiple R-squared: 0.0201, Adjusted R-squared: -0.08878
## F-statistic: 0.1846 on 1 and 9 DF, p-value: 0.6775Hydric Conclusion
There are basically no detectable relationships between water balance (CEWL, osml, or mass) and thermoregulation or microhabitat use in BNLL. The only significant relationship we found was between plasma osmolality and relative date of onset of estivation, but the effect size was miniscule.
Microhabitat Use ONLY
Prep Data
For logistic regression, I need to have the full, raw dataset (vs proportions by group for plots), so I will actually use the ‘microhab’ dataframe, created in the “Load Data” > “Microhbaitat Use” section.
Check Sex Differences
ggplot(MH_use_across_indivs_sex) +
geom_bar(aes(x = date_chunks,
y = MH_use_proportion,
fill = microhabitat
),
position = "fill",
stat = "identity") +
theme_classic() +
facet_wrap(~sex_M_F)
sex_micro_glm <- glm(data = microhab,
above_below ~ date_chunks*sex_M_F,
family="binomial")
summary(sex_micro_glm)##
## Call:
## glm(formula = above_below ~ date_chunks * sex_M_F, family = "binomial",
## data = microhab)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7145 -0.9537 -0.5474 0.9874 1.9861
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.33829 0.28088 -4.765 1.89e-06 ***
## date_chunksMay 10 - 20 1.05060 0.37044 2.836 0.00457 **
## date_chunksMay 21 - 31 0.50538 0.38811 1.302 0.19287
## date_chunksJun 1 - 11 0.78622 0.40228 1.954 0.05065 .
## date_chunksJun 12 - 22 1.88685 0.37360 5.050 4.41e-07 ***
## date_chunksJun 23 - Jul 3 2.48475 0.34300 7.244 4.35e-13 ***
## date_chunksJul 4 - 13 2.30337 0.40645 5.667 1.45e-08 ***
## sex_M_FMale -0.07878 0.39566 -0.199 0.84217
## date_chunksMay 10 - 20:sex_M_FMale -1.21595 0.54336 -2.238 0.02523 *
## date_chunksMay 21 - 31:sex_M_FMale -0.91084 0.54853 -1.661 0.09681 .
## date_chunksJun 1 - 11:sex_M_FMale -1.00676 0.54439 -1.849 0.06441 .
## date_chunksJun 12 - 22:sex_M_FMale -0.31563 0.48922 -0.645 0.51881
## date_chunksJun 23 - Jul 3:sex_M_FMale -0.60287 0.45270 -1.332 0.18295
## date_chunksJul 4 - 13:sex_M_FMale 0.32201 0.52066 0.618 0.53627
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2451.7 on 1768 degrees of freedom
## Residual deviance: 2060.7 on 1755 degrees of freedom
## AIC: 2088.7
##
## Number of Fisher Scoring iterations: 4drop1(sex_micro_glm)## Single term deletions
##
## Model:
## above_below ~ date_chunks * sex_M_F
## Df Deviance AIC
## <none> 2060.7 2088.7
## date_chunks:sex_M_F 6 2075.2 2091.2anova(sex_micro_glm)## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: above_below
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 1768 2451.7
## date_chunks 6 350.35 1762 2101.4
## sex_M_F 1 26.19 1761 2075.2
## date_chunks:sex_M_F 6 14.47 1755 2060.7probability_sex <- sex_micro_glm %>%
predict(microhab, type = "response")
predicted_binary_sex <- microhab %>%
mutate(pred_bin = if_else(probability_sex > 0.5, 1, 0))
ggplot(predicted_binary_sex) +
aes(x = date_chunks,
y = pred_bin,
color = sex_M_F) +
geom_jitter(alpha = 0.2) +
geom_smooth(method = "glm", method.args = list(family = "binomial"))## `geom_smooth()` using formula = 'y ~ x'
Date Only
ggplot(MH_use_across_indivs) +
geom_bar(aes(x = date_chunks,
y = MH_use_proportion,
fill = microhabitat
),
position = "fill",
stat = "identity") +
theme_classic()
micro_glm <- glm(data = microhab,
above_below ~ date_chunks,
family="binomial")
summary(micro_glm)##
## Call:
## glm(formula = above_below ~ date_chunks, family = "binomial",
## data = microhab)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6887 -0.8188 -0.6660 0.9292 1.7972
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.37845 0.19780 -6.969 3.19e-12 ***
## date_chunksMay 10 - 20 0.45773 0.26494 1.728 0.0841 .
## date_chunksMay 21 - 31 -0.01476 0.27166 -0.054 0.9567
## date_chunksJun 1 - 11 0.11857 0.26686 0.444 0.6568
## date_chunksJun 12 - 22 1.64082 0.23496 6.983 2.88e-12 ***
## date_chunksJun 23 - Jul 3 1.99483 0.21612 9.230 < 2e-16 ***
## date_chunksJul 4 - 13 2.52943 0.24567 10.296 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2451.7 on 1768 degrees of freedom
## Residual deviance: 2101.4 on 1762 degrees of freedom
## AIC: 2115.4
##
## Number of Fisher Scoring iterations: 4pairs(emmeans(micro_glm, "date_chunks"))## contrast estimate SE df z.ratio p.value
## (Apr 26 - May 6) - (May 10 - 20) -0.4577 0.265 Inf -1.728 0.5973
## (Apr 26 - May 6) - (May 21 - 31) 0.0148 0.272 Inf 0.054 1.0000
## (Apr 26 - May 6) - (Jun 1 - 11) -0.1186 0.267 Inf -0.444 0.9994
## (Apr 26 - May 6) - (Jun 12 - 22) -1.6408 0.235 Inf -6.983 <.0001
## (Apr 26 - May 6) - (Jun 23 - Jul 3) -1.9948 0.216 Inf -9.230 <.0001
## (Apr 26 - May 6) - (Jul 4 - 13) -2.5294 0.246 Inf -10.296 <.0001
## (May 10 - 20) - (May 21 - 31) 0.4725 0.256 Inf 1.843 0.5189
## (May 10 - 20) - (Jun 1 - 11) 0.3392 0.251 Inf 1.349 0.8284
## (May 10 - 20) - (Jun 12 - 22) -1.1831 0.217 Inf -5.448 <.0001
## (May 10 - 20) - (Jun 23 - Jul 3) -1.5371 0.197 Inf -7.818 <.0001
## (May 10 - 20) - (Jul 4 - 13) -2.0717 0.229 Inf -9.059 <.0001
## (May 21 - 31) - (Jun 1 - 11) -0.1333 0.258 Inf -0.516 0.9986
## (May 21 - 31) - (Jun 12 - 22) -1.6556 0.225 Inf -7.348 <.0001
## (May 21 - 31) - (Jun 23 - Jul 3) -2.0096 0.206 Inf -9.776 <.0001
## (May 21 - 31) - (Jul 4 - 13) -2.5442 0.236 Inf -10.760 <.0001
## (Jun 1 - 11) - (Jun 12 - 22) -1.5222 0.219 Inf -6.935 <.0001
## (Jun 1 - 11) - (Jun 23 - Jul 3) -1.8763 0.199 Inf -9.419 <.0001
## (Jun 1 - 11) - (Jul 4 - 13) -2.4109 0.231 Inf -10.440 <.0001
## (Jun 12 - 22) - (Jun 23 - Jul 3) -0.3540 0.154 Inf -2.301 0.2433
## (Jun 12 - 22) - (Jul 4 - 13) -0.8886 0.193 Inf -4.600 0.0001
## (Jun 23 - Jul 3) - (Jul 4 - 13) -0.5346 0.170 Inf -3.149 0.0273
##
## Results are given on the log odds ratio (not the response) scale.
## P value adjustment: tukey method for comparing a family of 7 estimatesanova(sex_micro_glm, micro_glm)## Analysis of Deviance Table
##
## Model 1: above_below ~ date_chunks * sex_M_F
## Model 2: above_below ~ date_chunks
## Resid. Df Resid. Dev Df Deviance
## 1 1755 2060.7
## 2 1762 2101.4 -7 -40.66probability_micro <- micro_glm %>%
predict(microhab, type = "response")
predicted_binary_micro <- microhab %>%
mutate(pred_bin = if_else(probability_micro > 0.5, 1, 0))
ggplot(predicted_binary_micro) +
aes(x = date_chunks,
y = pred_bin,
color = date_chunks) +
scale_color_manual(values = my_vir_cols_7) +
savs_ggplot_theme +
geom_jitter(alpha = 0.9) +
geom_smooth(method = "glm", method.args = list(family = "binomial"))## `geom_smooth()` using formula = 'y ~ x'
Pretty Barplot
par(mfrow=c(2,3))
ggplot(MH_use_across_indivs) +
geom_bar(aes(x = date_chunks,
y = MH_use_proportion,
fill = microhabitat
),
position = "fill",
stat = "identity") +
scale_fill_manual(values = my_brew_cols_4) +
savs_ggplot_theme
ggplot(MH_use_across_indivs) +
geom_bar(aes(x = date_chunks,
y = MH_use_proportion,
fill = microhabitat
),
position = "fill",
stat = "identity") +
scale_fill_manual(values = my_vir_cols_41) +
savs_ggplot_theme
ggplot(MH_use_across_indivs) +
geom_bar(aes(x = date_chunks,
y = MH_use_proportion,
fill = microhabitat
),
position = "fill",
stat = "identity") +
scale_fill_manual(values = my_vir_cols_42) +
savs_ggplot_theme
ggplot(MH_use_across_indivs) +
geom_bar(aes(x = date_chunks,
y = MH_use_proportion,
fill = microhabitat
),
position = "fill",
stat = "identity") +
scale_fill_manual(values = my_vir_cols_43) +
savs_ggplot_theme
ggplot(MH_use_across_indivs) +
geom_bar(aes(x = date_chunks,
y = MH_use_proportion,
fill = microhabitat
),
position = "fill",
stat = "identity") +
scale_fill_manual(values = my_vir_cols_44) +
savs_ggplot_theme
ggplot(MH_use_across_indivs) +
geom_bar(aes(x = date_chunks,
y = MH_use_proportion,
fill = microhabitat
),
position = "fill",
stat = "identity") +
scale_fill_manual(values = my_vir_cols_45) +
savs_ggplot_theme
votes: 1 2 3 4 5
~ temp interannual variation (Kat & Nicole ~ rainfall)
Notes:
- Kat tracked May - July & split micro use May/June/July ^ overlaid daytime and nighttime temps on the micro use fig
- Nicole tracked from May - mid-July & split micro use May/June/July
From Nicole’s paper: We then calculated the percent of time G. sila used each microhabitat in May, June, and July at each site. To compare the probability that a lizard would be found underground (in burrows) between the two sites, we ran a mixed-effects logistic regression model in R (R Core Team, 2020; RStudio, 2020, lme4 package v. 1.1-26, Bates et al., 2015) with time as a polynomial, site and month as fixed effects, and lizard ID as a random effect.